# This script calculates descriptive statistics and generates descriptive figures

# Last run 19 December 2024

dat_conjoint <- read_csv("data/dat_conjoint.csv")

# 1. Calculate survey respondent descriptives --------------

## 1.1 Country --------------

country_df <- dat_conjoint %>%
  select(response_id, country) %>%
  distinct() %>%
  mutate(N_variable = n()) %>%
  group_by(country) %>%
  summarize(N = n(),
            `N (\\%)` = paste0(n(), " (", round(n()/unique(N_variable)*100),"\\%)"),.groups = 'drop') %>%
  arrange(desc(N))

country_df <- bind_rows(country_df[1:9,],
                        tibble(country = c(paste0("Others (",nrow(country_df)-9,")"),"\\midrule Total"),
                               N = c(sum(country_df$N[-c(1:9)]),sum(country_df$N)),
                               `N (\\%)` =  c(paste0(sum(country_df$N[-c(1:9)]), " (", round(sum(country_df$N[-c(1:9)])/sum(country_df$N)*100),"\\%)"),sum(country_df$N)))) %>%
  select(-N)

addtorow <- list(); addtorow$pos <- list(nrow(country_df)); addtorow$command <- c("")
ptable(country_df,"output/country_df.tex")

## 1.2 Gender, location, position, years experience --------------
demographics_df <- dat_conjoint %>%
  select(response_id, position, location, years_experience, gender) %>%
  distinct() %>%
  select(-response_id) %>%
  pivot_longer(cols = everything(),names_to = "Variable",values_to = "Level") %>%
  mutate(Level = if_else(Level == "Academic & consultant","Academic \\& consultant", Level)) %>%
  group_by(Variable) %>%
  mutate(N_variable = n()) %>%
  group_by(Variable, Level) %>%
  summarize(`N (\\%)` = paste0(n(), " (", round(n()/unique(N_variable)*100),"\\%)"),.groups = 'drop') %>%
  mutate(Variable = str_replace_all(pattern = "_",replacement = " ",str_to_title(Variable)))

for(v in unique(demographics_df$Variable)){
  demographics_df$Variable[which(demographics_df$Variable==v)[-1]] <- ""
}
for(v in unique(demographics_df$Variable)){
  if(v %in% c("Gender","")==F){
    demographics_df$Variable[which(demographics_df$Variable==v)] <- paste0("\\midrule ",v)
  }
}

demographics_df$Level[which(is.na(demographics_df$Level))] <- "N/A"
demographics_df$Level <- paste0("\\multicolumn{1}{l}{",demographics_df$Level,"}")

addtorow <- list(); addtorow$pos <- list(nrow(demographics_df)); addtorow$command <- c("")
ptable(demographics_df,"output/demographics_1.tex")

## 1.3 Respondent carbon pricing views and country characteristics --------------

demographics_df <- dat_conjoint %>%
  select(response_id, cp_support, driver_domestic_international,driver_private_public,adaptation_mitigation,
         v2x_polyarchy, gdp_pc, oda_pc, cp_status_factor, tax_status_factor, ets_status_factor) %>%
  distinct() %>%
  pivot_longer(cols = -1,names_to = "Variable",values_to = "value") %>%
  group_by(Variable) %>%
  summarize(N = n(),
            Mean = round(mean(value,na.rm = T),digits = 2),
            `St. Dev` = round(sd(value,na.rm = T),digits = 2),
            Min = round(min(value,na.rm = T),digits = 2),
            Max = round(max(value,na.rm = T),digits = 2)) %>%
  mutate(Variable = case_when(Variable == "cp_support" ~ "Carbon pricing support",
                              Variable == "driver_domestic_international" ~ "International/domestic carbon pricing pressure",
                              Variable == "driver_private_public" ~ "Private/public carbon pricing demand",
                              Variable == "adaptation_mitigation" ~ "Adaptation/mitigation political priority",
                              Variable == "v2x_polyarchy" ~ "Electoral democracy",
                              Variable == "gdp_pc" ~ "Gross Domestic Product per capita (asinh)",
                              Variable == "oda_pc" ~ "Overseas Development Assistance per capita (asinh)",
                              Variable == "cp_status_factor" ~ "Carbon pricing status",
                              Variable == "tax_status_factor" ~ "Carbon tax status",
                              Variable == "ets_status_factor" ~ "ETS status")) %>%
  mutate(Variable = factor(Variable, levels = c("Carbon pricing support","International/domestic carbon pricing pressure","Private/public carbon pricing demand",
                                                "Adaptation/mitigation political priority","Electoral democracy","Gross Domestic Product per capita (asinh)",
                                                "Overseas Development Assistance per capita (asinh)","Carbon pricing status","Carbon tax status","ETS status"))) %>%
  arrange(Variable)

addtorow <- list(); addtorow$pos <- list(nrow(demographics_df)); addtorow$command <- c("")
ptable(demographics_df,"output/demographics_2.tex")

# 2. Carbon pricing policy adoption figures --------------

## 2.1 Prepare carbon pricing adoption/consideration data --------------
cpi <- read_excel("data/carbon_pricing_dashboard_data-latest.xlsx",skip = 1)
colnames(cpi)[c(3,4,5)] <- c("type","status","jurisdiction")
cpi <- cpi %>% 
  select(type, status, jurisdiction) %>%
  filter(is.na(type)==F)

### 2.1.1 Identify CPI country--------------
cpi$country <- ifelse(str_detect(cpi$type,pattern = "National"),cpi$jurisdiction,NA)
cpi$country[grep(cpi$jurisdiction,pattern = "EU")] <- "European Union"
cpi$country[grep(cpi$jurisdiction,pattern = "RGGI")] <- "United States of America"

# Geocode jurisdiction name
cpi_geocoded <- cpi %>%
  filter(is.na(country)) %>%
  geocode(jurisdiction,method = "osm",lat = latitude, long = longitude) %>%
  reverse_geocode(lat = latitude, long = longitude, method = "osm", address = address_found, full_results = TRUE)

cpi_geocoded <- cpi_geocoded %>%
  select(jurisdiction, country_code) %>%
  mutate(country = countrycode(country_code,origin = "iso2c",destination = "country.name"))

cpi_geocoded$country[which(cpi_geocoded$jurisdiction=="Guangdong (except Shenzhen)")] <- "China"

cpi$country[which(is.na(cpi$country))] <- cpi_geocoded$country
cpi$country[which(cpi$jurisdiction=="Taiwan, China")] <- "Taiwan"

# Regularize country names
cpi$country[which(cpi$country!="European Union")] <- countrycode(cpi$country[which(cpi$country!="European Union")],origin = "country.name",destination = "country.name")

# Mark advanced country CPIs
advanced <- c("Australia","Austria","Canada","Denmark","Estonia","European Union",
              "Finland","France","Germany","Hungary","Iceland","Ireland","Israel","Japan","Latvia","Liechtenstein",
              "Luxembourg","Netherlands","New Zealand","Norway","Poland","Portugal","Russia","Slovenia","Spain",
              "South Korea","Sweden","Switzerland","Taiwan","United Kingdom","United States")

cpi$developing <- if_else(cpi$country %in% advanced,0,1)

cpi$adopted_year <- as.integer(sapply(str_extract_all(cpi$status,pattern = "[0-9]{4}"),function(x){x[1]}))
cpi$abolished_year <- as.integer(sapply(str_extract_all(cpi$status,pattern = "[0-9]{4}"),function(x){ifelse(length(x)==2,x[2],NA)}))
cpi$status <- ifelse(str_detect(cpi$status, pattern = "Implemented in [0-9]{4}$"),"Implemented",
                     ifelse(str_detect(cpi$status, pattern = "Under development"),"Implemented",
                            ifelse(str_detect(cpi$status, pattern = "Under consideration"),"Under consideration",
                                   "Abolished")))

## 2.2 CPI adoptions-by-year plot --------------
# Calculate number of implemented instruments in each year, by developing/advanced
cpi_year <- cpi %>%
  filter(status == "Implemented") %>%
  group_by(developing,adopted_year) %>%
  summarize(n=n(),.groups = 'drop') %>%
  group_by(developing) %>%
  mutate(n = cumsum(n)) %>%
  rename(year = adopted_year)

idx <- which(cpi$status == "Abolished")

# Add in abolished instruments during years of implementation
for(i in idx){
  year_idx <- which(cpi_year$developing==cpi$developing[i] & cpi_year$year %in% cpi$adopted_year[i]:(cpi$abolished_year[i]))
  cpi_year$n[year_idx] <- cpi_year$n[year_idx] + 1
}

# Calculate number of instruments under consideration
cpi_year <- bind_rows(cpi_year,
                      tibble(developing = c(0,1),
                             year = max(cpi_year$year) + 2,
                             n = c(length(which(cpi$developing==0 & cpi$status=="Under consideration")),
                                   length(which(cpi$developing==1 & cpi$status=="Under consideration")))))

# Interpolate years
cpi_year <- cpi_year %>% filter(year >= 2000)

cpi_year_plotdat <- tibble(developing = NA_integer_, year = NA_integer_, n = NA_integer_, .rows = 0)
  
for(yr in (2000:2026)){
  for(dv in c(0,1)){
    if(yr %in% cpi_year$year[which(cpi_year$developing==dv)] == F){
      if(yr >= min(cpi_year$year[which(cpi_year$developing==dv)])){
        cpi_year_plotdat <- bind_rows(cpi_year_plotdat,
                              tibble(developing = dv,
                                     year = yr,
                                     n = cpi_year$n[which(cpi_year$year==max(cpi_year$year[which(cpi_year$developing==dv & cpi_year$year<=yr)]) & cpi_year$developing==dv)]))
      }
    } else {
      cpi_year_plotdat <- bind_rows(cpi_year_plotdat,
                                    cpi_year[which(cpi_year$year==yr & cpi_year$developing==dv),])
    }
  }
}

cpi_year_plotdat <- cpi_year_plotdat %>% 
  filter(year!=2025) %>%
  arrange(developing,year)

# Create country detail label
country_detail_label <- paste0(names(table(cpi %>% filter(status == "Implemented" & developing == 1) %>% select(country))),collapse = ", ")
country_detail_label <- paste0("<span style='color:#1E2F57;'>**Developing countries with carbon pricing implemented or scheduled:**</span><br>", country_detail_label)
country_detail_label <- gsub(pattern = ", Montenegro",replacement = ",<br>Montenegro",country_detail_label)

cpi_year_plot <- ggplot(cpi_year_plotdat, aes(x = year, y = n, group = developing)) +
  geom_col(color = "black", width = 0.75, fill = c(rep("#883369",26),rep("#1E2F57",18)),alpha = c(rep(1,25),0.5,rep(1,17),0.5)) +
  labs(x = NULL, y = "Number of regional, national, and subnational carbon pricing initiatives") +
  scale_y_continuous(expand = c(0,0), limits = c(0,105), breaks = seq(0,100,by=20)) +
  scale_x_continuous(expand = expansion(add = c(0,1)),breaks = 2000:2026,labels = c(2000:2024,"","Under<br>Consideration"))+
  annotate("richtext",x = c(2024.25,2024.25), y = c(93,38),label = c("**Developed**","**Developing**"),hjust = 1,
           color = c("#883369","white"), label.color = c("#883369","white"), linewidth = 1, fill = c("white","#1E2F57"),
           size = 4.5) +
  annotate("shadowtext",x = c(2024,2024), y = c(88,32),label = c("56","34"),hjust = 0.5,
           color = "white",bg.color = "black",
           size = 4.5) +
  annotate("shadowtext",x = c(2026,2026), y = c(28,19),label = c("9","21"),hjust = 0.5,
           color = "white",bg.color = "black",
           size = 4.5) +
  annotate("richtext",x = 2001, y = 90,hjust = 0,label.color = "#1E2F57",label.size = 0,
           label = country_detail_label, size = 4.5) +
  theme_classic() +
  theme(axis.text = element_markdown(face = c(rep("plain",25),"bold"),size = 10),axis.ticks.x = element_blank(),
        axis.title.y = element_text(face = "bold",size = 12))

ggsave(cpi_year_plot,file = "output/cpi_year_plot.png",width = 12,height = 8,units = "in",dpi = 320)

## 2.3 CPI adoption/consideration map --------------

# Create dataframe of developing countries and maximum CPI pricing status (max = adopted anywhere in country)
map_plotdata <- cpi %>%
  mutate(type = case_when(status=="Under consideration" & developing == 1 ~ "Under consideration",
                          status=="Under consideration" & developing == 0 ~ "Under consideration_developed",
                          (status=="Implemented" | status=="Scheduled") & developing == 1 ~ "Adopted",
                          (status=="Implemented" | status=="Scheduled") & developing == 0 ~ "Adopted_developed"),
         type = factor(type, levels = c("Under consideration","Adopted","Under consideration_developed","Adopted_developed"))) %>%
  select(country, type) %>%
  na.omit() %>%
  group_by(country) %>%
  mutate(keep = if_else(as.numeric(type) == max(as.numeric(type)),1,0)) %>%
  filter(keep == 1) %>%
  select(-keep) %>%
  distinct() %>%
  # Drop developed countries considering carbon pricing
  filter(type != "Under consideration_developed")

# Replace European Union with constituent countries
eu <- c("Iceland", "Liechtenstein", "Norway", "Austria", "Belgium", "Bulgaria", "Croatia", "Cyprus", "Czech Republic", 
        "Denmark", "Estonia", "Finland", "France", "Germany", "Greece", "Hungary", "Ireland", "Italy", "Latvia", 
        "Lithuania", "Luxembourg", "Malta", "Netherlands", "Poland", "Portugal", "Romania", "Slovakia", "Slovenia", "Spain", "Sweden")

map_plotdata <- bind_rows(map_plotdata,
                          tibble(country = eu[which(eu %in% map_plotdata$country==F)], type = "Adopted_developed"))

map_plotdata <- map_plotdata %>%
  rename(region = country) %>%
  mutate(region = if_else(region == "Bosnia & Herzegovina","Bosnia and Herzegovina",
                          if_else(region == "Côte d’Ivoire","Ivory Coast",
                                  if_else(region == "United Kingdom","UK",
                                          if_else(region == "United States","USA",region)))))

world_map <- map_data("world")

map_plotdata <- map_plotdata %>%
  bind_rows(world_map %>%
              select(region) %>%
              filter(region %in% map_plotdata$region==F) %>%
              distinct() %>%
              mutate(type = "Developed")) %>%
  mutate(type = if_else(region=="Antarctica",NA,type))

map_plot <- ggplot(map_plotdata, aes(map_id = region)) +
  geom_map(aes(fill = type), color = "grey30",map = world_map,linewidth = 0.25) + 
  expand_limits(x = world_map$long, y = world_map$lat) +
  coord_sf(default_crs = 4326,
           ylim = c(-50, 90)
  ) + 
  scale_fill_manual(values = c("#33A02C","grey30","grey90","#B2DF8A")) +
  theme_map() + 
  theme(legend.position = "None")

ggsave(map_plot,file = "output/map_plot.png",width = 12.36,height = 10.02,units = "in")

